home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / linalg / multiply.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  3.4 KB  |  138 lines
  1. /* linalg/multiply.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_linalg.h>
  26.  
  27. #define SWAP_SIZE_T(a, b)  do { size_t swap_tmp = a; a = b; b = swap_tmp; } while(0)
  28.  
  29. int
  30. gsl_linalg_matmult (const gsl_matrix * A, const gsl_matrix * B, gsl_matrix * C)
  31. {
  32.   if (A->size2 != B->size1 || A->size1 != C->size1 || B->size2 != C->size2)
  33.     {
  34.       GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN);
  35.     }
  36.   else
  37.     {
  38.       double a, b;
  39.       double temp;
  40.       size_t i, j, k;
  41.  
  42.       for (i = 0; i < C->size1; i++)
  43.     {
  44.       for (j = 0; j < C->size2; j++)
  45.         {
  46.           a = gsl_matrix_get (A, i, 0);
  47.           b = gsl_matrix_get (B, 0, j);
  48.           temp = a * b;
  49.           for (k = 1; k < A->size2; k++)
  50.         {
  51.           a = gsl_matrix_get (A, i, k);
  52.           b = gsl_matrix_get (B, k, j);
  53.           temp += a * b;
  54.         }
  55.           gsl_matrix_set (C, i, j, temp);
  56.         }
  57.     }
  58.  
  59.       return GSL_SUCCESS;
  60.     }
  61. }
  62.  
  63.  
  64. int
  65. gsl_linalg_matmult_mod (const gsl_matrix * A, gsl_linalg_matrix_mod_t modA,
  66.             const gsl_matrix * B, gsl_linalg_matrix_mod_t modB,
  67.             gsl_matrix * C)
  68. {
  69.   if (modA == GSL_LINALG_MOD_NONE && modB == GSL_LINALG_MOD_NONE)
  70.     {
  71.       return gsl_linalg_matmult (A, B, C);
  72.     }
  73.   else
  74.     {
  75.       size_t dim1_A = A->size1;
  76.       size_t dim2_A = A->size2;
  77.       size_t dim1_B = B->size1;
  78.       size_t dim2_B = B->size2;
  79.       size_t dim1_C = C->size1;
  80.       size_t dim2_C = C->size2;
  81.  
  82.       if (modA & GSL_LINALG_MOD_TRANSPOSE)
  83.     SWAP_SIZE_T (dim1_A, dim2_A);
  84.       if (modB & GSL_LINALG_MOD_TRANSPOSE)
  85.     SWAP_SIZE_T (dim1_B, dim2_B);
  86.  
  87.       if (dim2_A != dim1_B || dim1_A != dim1_C || dim2_B != dim2_C)
  88.     {
  89.       GSL_ERROR ("matrix sizes are not conformant", GSL_EBADLEN);
  90.     }
  91.       else
  92.     {
  93.       double a, b;
  94.       double temp;
  95.       size_t i, j, k;
  96.       size_t a1, a2, b1, b2;
  97.  
  98.       for (i = 0; i < dim1_C; i++)
  99.         {
  100.           for (j = 0; j < dim2_C; j++)
  101.         {
  102.           a1 = i;
  103.           a2 = 0;
  104.           b1 = 0;
  105.           b2 = j;
  106.           if (modA & GSL_LINALG_MOD_TRANSPOSE)
  107.             SWAP_SIZE_T (a1, a2);
  108.           if (modB & GSL_LINALG_MOD_TRANSPOSE)
  109.             SWAP_SIZE_T (b1, b2);
  110.  
  111.           a = gsl_matrix_get (A, a1, a2);
  112.           b = gsl_matrix_get (B, b1, b2);
  113.           temp = a * b;
  114.  
  115.           for (k = 1; k < dim2_A; k++)
  116.             {
  117.               a1 = i;
  118.               a2 = k;
  119.               b1 = k;
  120.               b2 = j;
  121.               if (modA & GSL_LINALG_MOD_TRANSPOSE)
  122.             SWAP_SIZE_T (a1, a2);
  123.               if (modB & GSL_LINALG_MOD_TRANSPOSE)
  124.             SWAP_SIZE_T (b1, b2);
  125.               a = gsl_matrix_get (A, a1, a2);
  126.               b = gsl_matrix_get (B, b1, b2);
  127.               temp += a * b;
  128.             }
  129.  
  130.           gsl_matrix_set (C, i, j, temp);
  131.         }
  132.         }
  133.  
  134.       return GSL_SUCCESS;
  135.     }
  136.     }
  137. }
  138.